One-dimensional counterion gas between charged surfaces: 
Exact results compared with weak- and strong-coupling analysis 
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We evaluate exetctly the statistical integral for an inhomogeneous one- dimensional counterion-only 
Coulomb gas between two charged boundaries and from this compute the effective interaction, or 
disjoining pressure, between the bounding surfaces. Our exact results are compared with the limit- 
ing cases of weak and strong coupling which are the same for ID and 3D systems. For systems with 
a large number of counterions it is found that the weak coupling (mean-field) approximation for 
, the disjoining pressure works perfectly and that fluctuations around the mean-field in ID are much 

O ' smaller than in 3D. In the case of few counterions it works less well and strong coupling approx- 

imation performs much better as it takes into account properly the discreteness of the counterion 
charges. 

^ \ PACS numbers: 05.20.Jj, 61.20.' 

+-J , Keywords: 

I. INTRODUCTION 

Electrostatic effects are predominant in a variety of soft condensed matter systems such as polymers in solutions, 
colloidal suspensions and much of the physics of membranes and films [H, The statistical mechanics of these 
problems cannot in general be solved explicitly; however a number of approaches have been devised, starting with the 
mean- field Poisson-Boltzmann equation, that allow for an approximate solution to the problem Q . One can improve 
the accuracy of mean-field calculations by looking at the fluctuations about the Poisson-Boltzmann solution, and in 
a number of cases these fluctuation effects are essential to capture the quantitative behavior of the systems under 
study. The approach of this type can be shown to be valid in the so called weak-coupling regime. More recently the 
strong-coupling expansion, in its essence similar to the virial expansion, has been developed and successfully applied 
to a number of interesting problems [!, 13] . Despite the success of these approaches many questions remain to be 
answered about the intermediate regime and how the cross-over between the weak- and strong-coupling limits occurs. 
Here we show that the one-component ID inhomogeneous Coulomb fluid can be analytically solved and thus provides 
an ideal test-bed to study the domains of validity of these commonly used approximation schemes. 

For a system of two charged surfaces interacting in an electrolyte solution, the effective coupling strength may 
depend on the inter-plate separation, and so the calculation of the force between the plates requires an understanding 
of the both weak and strong coupling regimes and the transition between them. In two important papers by Lenard 
Q and Edwards and Lenard Q a complete solution of the thermodynamics of the two component Coulomb gas in 
one dimension was given. The second of these papers was based on the fact that the underlying Sine-Gordon theory 
for this system can be mapped onto quantum mechanics in one dimension and the thermodynamics can be solved in 
terms of the ground state energy of the corresponding Schrodinger equation. Subsequently the system was studied in 
the presence of an electric field and it was shown that it always behaves as a dielectric because of a dimerization of 
the individual charge carriers Q- When the system is not in an electroneutral state, i.e. when the surface charges 
corresponding to the applied field are not exactly compensated by the counterions of the system, the physics of 
the system is subtly different and confinement phenomena between oppositely charged particles appear [8|, l9( . This 
confinement is directly related to the fact that the system is dielectric and non-conducting. More recently, the one- 
dimensional Coulomb gas thermodynamics has been studied taking into account additional hard-core and dipolar 
interactions (To| . The model has also been studied in the finite size setting of a soap film type model where charge 
regulation mechanism model of the surface, in the presence of salt solution (ions and counterions), is included [llj . 

In the present paper we study a related model to those outlined above: a one-dimensional model of counterions 
confined between two plates, carrying fixed but not necessarily equal charges. This situation is quite different to most 
of those studied above as most of them have concentrated on symmetric electrolyte systems. We place particular 
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emphasis on the computation of the pressure of the system which is the effective interaction between the two plates 
modified by the presence of counterions. The effect of asymmetric charge distribution on the plates is analyzed 
in detail. This system, while being somewhat idealized, is a simple model of an array of charged smectic bilayers 
sandwiched between two fixed parallel charged plates. If the charge on the bilayers is uniform then the use of the 
one-dimensional Coulomb interaction is justified. 

Another reason for studying this simple system is to check the validity of the various approximation schemes 
alluded to above. The mean-field and strong-coupling approximations in the planar geometry are independent of 
the dimensionality of the problem. This is due simply to the fact the both Poisson-Boltzmann theory as well as 
the strong-coupling theory are one-dimensional effective theories since the corresponding fields depend only on the 
coordinate parallel to the bounding surfaces normals. However, fluctuations about the mean field, or the Poisson- 
Boltmzann configuration, do depend on the dimensionality of the problem as they correspond to thermal Casimir 
or zero-frequency van-der-Waals forces. This makes these model calculations particularly appealing since they can 
clarify the role the fluctuations play in Coulomb systems. 

The paper is organized as follows: first we use the Edwards-Lenard path integral formulation for the problem 
adapted to apply to a finite electroneutral system with surface charges. The exact expression obtained for the force 
between the two plates is evaluated numerically and compared with the results of Monte-Carlo simulation for a wide 
range of cases: symmetric/asymmetric surface charges, number of counterions, etc. We find excellent agreement 
between theory and simulation in all cases as we should expect since the theoretical solution is exact. In the next 
section we compare the results with the Poisson-Boltzmann mean-field approximation and investigate its domain of 
validity. We present two methods for calculating the effects of field fluctuations about the mean-field solution and 
compare with our exact results. The strong-coupling approximation is then applied to the system and a similar 
comparison with the exact results is made. 



II. THE COUNTERION MODEL 



We consider a system of charged particles interacting via a Coulomb potential in one dimension. It can also be 
considered as a system of uniformly charged infinite two dimensional sheets in three dimensions, where the sheets are 
perpendicular to the direction x, arc uniformly charged and can only move in this direction. The Hamiltonian for 
this system is given by 

0n = -—22z a Z/3\x a -X/3\, (1) 

a,0 

where, e is the unit electric charge, /? = 1/ksT, z a is the valence of particle (or sheet) a and x a is its position. We 
will consider systems in a canonical formulation which are overall charge neutral, so that 

5>a = 0" (2) 

a 

It is convenient to rewrite the Hamiltonian as 

/3e 2 \ - 

(3H = —22z a z p [x a +x fj - 2min(x Q ,a; / 3)] 

= ^Y^2 z a z /3 rmn(x a ,xp), (3) 

where we have used charge neutrality to go from the first to second equation above. We may now write the Hamiltonian 
as an expectation value over a "Brownian motion" 



where ip is Brownian motion with correlation function 



^z a il>(x a ) 



(4) 



(ip(x)ip(y)) = /3e 2 min(x,y), 



(5) 
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and started at the value zero at x = 0, so 

i>(0) = 0. (6) 
The Boltzmann weight for any configuration can thus be written as 

exp(-/3H) = (exp[i}^z a j)(x a )]). (7) 

a 

In path-integral notation the measure on the Brownian motion ip as defined by Eq. is given by 



</>(£) 



(0[i>]) = / / c#]exp -— - / dx -^J. 0[ip], 



•0(0) 



2e 2 P J V dx 



dip(x) 



(8) 



where L is the overall length of the system. 

The sum over a has three distinct contributions: that from the counterions with 1 < a < N, for which from now 
on we assume have valence z a = 1, and those from the surface charge at x = 0, with valence z = — (N + P)/2, 
and from the surface charge at x^v+i = L, with valence z^+i = — (N — P)/2. The system is globally electroneutral, 
consisting of two surface charges and their counterions. For the moment, we will assume that N ± P are even integers, 
which will be the case if the counterions originate by being released into solution from initially neutral surfaces; i.e., 
each surface charge is an integer multiple of the unit of charge of the counterions. The generalization of the approach 
due to relaxing the latter assumption is given in the next section. We then have 

A N + P , N-P 



£ ZcMXa) = ]T — ±— m — ^(L), (9) 

a i—1 

and the Boltzmann weight for a given configuration is thus 

BW({x t }) = ^(^^(x0-^m^^-^(L)^^^ ^. (10) 
The partition function is then obtained as 

Zn = £ IidxiBW{{ Xi }^ 

(L \ N 
jf dxexp(ii>(x))j exp(-^(0)^±^-#(L)t^jy (11) 

We now introduce an arbitrary fugacity k having the dimensions of inverse length and write 
k n Z n = / eX p(y#(0)^^-^(Z)^^ 

dx exp(iip(x)) 



iifj(L) — iXNjcxp^Kj dxexp(i\ + itp(x)) j J, (12) 

although, clearly, the value of k should not affect the final physical results since charge neutrality will enforce the 
constraint N' — N. This is ensured by the integral over dX. If we define a shifted Brownian motion, 

4>(x) = A + tp(x), (13) 

the starting position of the field (f> at x = is simply (f)(0) = A but the integration over the end point remains free 
and we have 



J N 



'"^^exp ( -i(/>(0)E±f. - i(j) (L)^LJ^) exp [k [ dv cxpl io[ v ) ) \ ). (14) 
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The electroneutrality constraint has thus been absorbed into the integration over the starting point of the shifted field 

It is convenient, for what comes later, to generalize the result to the case where the counterions have valence 
z a = q, 1 < a < N . This is easily accomplished by replacing e by qe in the measure for the Brownian motion in Eq. 
([8]). This follows from the forgoing analysis and rescaling the field c/>(x) — > <j>(x)/q. In path integral notation we then 
have 



Z N — 



2 71 



#(0) 
2tt 



d(f>(L) exp -i</>(0) 



N + P 



i<t>{L) 



N - P 



0(0) 



d[4>] exp 



dx 



1 fd4{x 
2q 2 e 2 f3\ dx 



k exp (i(f>(x)) 



(15) 



as our final functional integral expression for the partition function. 



III. EXACT EVALUATION OF THE PARTITION FUNCTION 



We use the Feynman formula, which is the Euclidean version of the equivalence between path integrals and the 
propagator for the Schrodinger equation. We have that 



do f<t>(L) I pL 

dcj){L) J d[(j>] exp I — / dx 

CO J 7. \ JO 



1 



2q 2 e 2 /3 \dx 



— k exp (icf>(x)) 



exp I —i(f)(L) 



N -P 



e~xp(—LH) exp 



N -P 



where q is the counterion valence and where H is the complex operator/Hamiltonian given by 

q 2 e 2 (3 d 2 



H 



The partition function can thus be expressed as 



2 dz 2 



k exp(zz) 



kTZ n {L) 



27r dz ( . N + P \ « ( . N-P 

— exp y-iz — - — I exp(-LH) exp y~iz — - — 



(16) 



(17) 



(18) 



We should note that the operator H is complex in our problem but in the case of symmetric electrolytes [6j it is 
real. Because of this we analyse the problem via Fourier analysis rather than on the basis of the eigenfunctions of the 
operator H. 

We define integers Mi — (N + P)/2 and A/2 = (N — P)/2 and the surface charges are —M\qe and —M^qe at x = 
and x = L, respectively, with M± + M2 — N . A quick preliminary check of our approach is to consider the perfect 
gas limit where e = 0. In this case we find trivially that 



k n Zn(L) = K 



N 



£ / Mi+M2 

(Mi +M 2 )! 



and thus 



L N 

Zn(L) = — , 

which is the perfect gas result. In order to proceed further we consider the evaluation of 

f(z;L)=exp(-LH)f(z;Q). 



(19) 



(20) 



(21) 



We can write [11 1 



f(z) — a(n, L) exp(inz) 



(22) 
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and the induced evolution equation for the Fourier coefficients is 

da(n,L) n 2 q 2 e 2 f3 



dL 



a(n,L) + na(n - l,L) . (23) 



In our problem we have f(z, M2;0) = exp (— izM 2 ) and thus the initial condition on the Fourier coefficients 
a(n, M 2 ; 0) = 8 n -M 2 - We define a n = n n b n and find that 

db(n,M 2 ;L) n 2 q 2 e 2 (3, / r . ,, rX 

V rfL 7 = \^ Kn, M 2 ; L) + 6(n - f , Af 2 ; L) , (24) 

with the initial condition b(n, M 2 ; L) = n n 5„^M 2 an d from Eq. (|18|) we find that we can express the partition function 
for system of unit charges M% and M 2 on the left and right boundaries respectively as 

Z Ml)Ma (L) = b(M 1 ,M 2 ;L) 1 (25) 

but with the initial condition b n (0, M 2 ) = 5 n .-M 2 \ as we expect, the fugacity k, which was introduced on dimensional 
grounds, does not enter in the final physical result. The evolution Eq. (f24|) can be written as 

1 db(M u M 2 ;L) M 2 q 2 e 2 Z Ml -i,M a , 2fi s 

Pm u m 2 pb{MijM2 . L) dL 2 + f3z MuM2 ' 1 bj 

where Pm-l,M 2 is the pressure of the system. Clearly, the last term makes a positive contribution to the pressure as 
it is the ratio of two partition functions for systems with different surface charges. Equation (|2"rJ)) has a nice physical 
interpretation. The second term is, up to the factor of /3, the average density of counterions at z = because it is the 
restricted partition function where at least one particle is on the surface z = (thus giving a reduction in the surface 
charge to M — 1) normalized by the full partition function. We can thus write 

PMuM* = -y + 4p(0) , (27) 

where p{x) — (X^i ${ x ~ x i)) i s the average value of the counterion density at the surface z = and o\ — —Miqe 
the corresponding surface charge. However, this is simply the contact-value theorem for electrostatic systems known 
to be exact in any dimension [lO . [l2l . [l3l Il4| . In fact the contact- value theorem can be demonstrated via an extension 
of the path integral methods used here to higher dimensional systems with (parallel) planar geometries [151 ] . The 
average counterion density at the point x can be shown [(| [ll[ to be equal to 



p{x) = (exp(icf)(x))) 

f 2 * dz 

/ — exp (— iM\z) evp(-xH) exp(iz) exp(— (L — x)H) exp (— iM 2 z) . 
Jo 2lT 



_ __1 f 2 * d 

Zm x ,m 2 Jo 

If we set x = in this formula we obtain 



P(0) = , (28) 



in agreement with the previous discussion. 



A. The general formalism 



The charges on the surfaces at z = and L can be denoted <j\ and a 2 respectively. Without loss of generality we 
may assume that q > and \a 2 \ > |cti|, which allows a\ and a 2 to be of opposite sign, and we can define the charge 
asymmetry parameter £ = u 2 j<j\ with — 1 < £ < 1. All other cases can be mapped onto this interval with appropriate 
rescaling of the parameters. Note that the values £ = 1 and £ = — 1 represent special cases of symmetric system with 
<j\ = cr 2 , and antisymmetric system with a± = — <j 2 and no counterions between surfaces (which thus reduces to the 
trivial case of a planar capacitor). 

For the above analysis we see that 

N - P 

C = , (29) 

S N+P ' V ; 
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<7i = —M%qe and a 2 — —M2qe. The independent parameters are conveniently chosen as a\,C,N. With M x and M2 
integers, the analysis so far allows only for particular discrete values of C- However, the approach can be extended 
to accommodate all values of ( and we state the generalization of the algorithm here. Charge neutrality determines 
that the valence q of the counterfoil satisfies 



iVe 



(30) 



We are allowing that q need not be an integer. 

In addition to the Fourier coefficients b(n, M2; L) we introduce Fourier coefficients c(n, Mi; L). Equation (|24|) gives 
the evolution of the b(n, M2; L) from z = 0, the surface with assigned charge ai. The coefficients c(n, Mi; L) obey a 
similar evolution equation but now with the surface charge assignment reversed; the charge at z = is now 02- These 
two sets of coefficients correspond to complementary approaches; the set {&} describes evolution from the surface with 
charge ci to that with charge a 2 , and vice- versa for the set {c}. 

We define 



f 



1 + C 

The evolution equations now become 

db(n,M 2 ;L) _ (n - i) 2 f 

dL ~ 2 
dc(n,Mi;L) _ {n-rnf 



M x = Int(aiV) , r]i = aN -Mi, rj 2 = 1 - rji , M 2 = N -M x -1 



dL 



pq 2 e 2 b(n, M 2 ; L) + b(n - 1, M 2 ; L) , 
Pq 2 e 2 c(n, M X ;L) + c(n - I, M X ;L) . 



(31) 

(32) 
(33) 



The initial conditions are 6(n, M 2 ; L = 0) = 8 n , c(n, Mi; L = 0) = 5 n _Mi ■ 

Then the partition function is given by alternative formulas Z(<ti, (, N; L) = b(M\ + I, M 2 ; L) — c(M2 + 1, Mi; L), 
and the pressure is given equivalently by 



r {Mi + 1 - m ? 



2 2 b(M u M 2 ;L) 
q e + 



P(L) = { 



(M 2 + i- m ) 2 



2 2 

q e 



Z(ai,(,N;L) 

c(M 2 ,Mi;L) 
Z(ai,t,N;L) 



(34) 



The partition function Z(ai, C, N; L) is also given by 



M 2 + l 



Z(ai,C,N;L)= ^ b(-j + 1, M 2 ; L - x) c{j,Mi;x), 
i=-Mi 



(35) 



for any x, < x < L. This acts as a check on the numerics. 
The counterion number density p{x; ai, C, TV, L) is given by 



p(x;ai,C,N, L) 



1 



m 2 +i 



Z(ai,(,N;L) 



J2 b(-j,M 2 ;L-x)c(j,M_ 



i;x) 



j=-Mi 



and another check on the numerics is that 



dxp{x;ai,C, 1 N,L) = N, 



(36) 



(37) 



i.e. the system is electroneutral. 



IV. APPROXIMATE EVALUATIONS OF THE PARTITION FUNCTION 



A traditional approach to the study of charged (bio)colloidal systems is the mean-field Poisson-Boltzmann (PB) 
formalism which is applicable for weak surface charges, low counter-ion valency and high temperature [1 61 ] - The limi- 
tations of this approach are evident when applied to highly-charged systems where counterion-mediated interactions 
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between charged bodies start to deviate substantially from the accepted mean-field wisdom 0, 0| ■ One of the recent 
fundamental advances in this field has been the systematization of these non-PB effects based on the notions of weak 
and strong coupling approximations. The latter approach has been pioneered by Rouzina and Bloomfield [13, elabo- 
rated later by Shklovskii et al. [H[, Levin et al. [ljjj], and brought into final form by Netz et al. [U, 0, H3, [Hf- These 
two approximations allow for an explicit and exact treatment of charged systems at two disjoint limiting conditions 
whereas the parameter space in-between can be analyzed onl y appro ximatel y I20L [2TL I22I I23I HH [25[ and is mostly 
accessible solely via computer simulations @, H H3, E], [H, H Hf HI HI HI IMllt 

Both the weak- and the strong-coupling approximations are based on a functional integral or field-theoretic rep- 
resentation [H, [33| of the grand canonical partition function for a system composed of fixed surface c harg es with 
intervening mobile counterions, and depend on the value of a single dimensionless coupling parameter S [2Ct l2l|. In 
three dimensions, S is proportional to the ratio between two relevant length scales, namely, the Bjerrum length and 
the Gouy-Chapman length. The Bjerrum length, defined as £b = e 2 /(47T££o/cbT), is the distance at which two unit 
charges interact with thermal energy k^T (in water at room temperature, one has £b — 0.7nm). If the charge valency 
of the counterions is q then the corresponding length scales as q 2 £B- Similarly, the Gouy-Chapman length, defined as 
fi — e/ (27rg^B|c|), is the distance at which a counterion interacts with a macromolecular surface (of surface charge 
density a) with an energy equal to k^T. The 3D electrostatic coupling parameter measures the competition between 
ion-ion and ion-surface interactions and is given by [20L l2~l| 



S = S 3D = q 2 £ B /fi = 27rg 3 4|a|/e . (38) 

Physically, the weak-coupling (WC) regime 5 <C 1 (appropriate for low valency counterions and/or weakly charged 
surfaces), is characterized by the fact that the width, fj,, of the counterion layer near the surfaces is much larger than 
the separation between two neighboring counterions in solution, and thus the counterion layer behaves basically as a 
three-dimensional gas. Each counterion in this case interacts with many others and the collective mean-field approach 
of the Poisson-Boltzmann (PB) type is completely justified. 

On the other hand in the strong-coupling (SC) regime H > 1 (appropriate for high valency counterions and/or 
highly charged surfaces), the mean distance between counterions, a± — \o~\/qe, is much larger than the layer width 
(i.e., a±/ fx ~ v3 3> 1), indicating that the counterions are highly localized laterally and form a strongly correlated 
quasi-two-dimensional layer next to a charged surface. In this case, the weak-coupling approach breaks down due to 
strong counterion-surface and counterion-counterion correlations. Since counterions can move almost independently 
from the others along the direction perpendicular to the surface, the collective many-body effects that enable a 
mean-field description are absent, necessitating a complementary SC description [20l. [2l|. 

Formally, the WC limit can be identified with the saddle-point approximation of the field theoretic representation 
of the grand canonical partition function, and reduces to the mean-field PB theory at the lowest order for H — > 0. The 
quadratic fluctuations around the mean field provide a second-order correction to the mean-field solution for small 
H < 1 [13, [H, H [H, [H, US H, HI. The SC approximation has no PB-like collective mean-field [U Hl[ since it is 
formally equivalent to a single particle description obtained from a systematic 1/5 expansion in the limit 5 — * 00, 
and corresponds to two lowest order terms in the virial expansion of the grand canonical partition function. The 
consequences and the formalism of these two limits of the Coulomb fluid description have been explored widely and 
in detail (for reviews, see Refs. 0,[3|)- 

The concept of weak- and strong-coupling electrostatics can be easily generalized to other dimensions. In one 
dimension (ID), the Coulomb interaction between two unit charges may be written as vid(x) = — \x\ksT / £b, where 
£b = 2fcsT/e 2 may be regarded as a ID Bjerrum length. Likewise, the interaction strength with the boundary charges 
in an asymmetric system may be characterized by the ID Gouy-Chapman lengths /ii = £Be/q\a\\ and /12 = £ B e/q\a2\- 
We may proceed by using only ji = fix = £b&i '<?|ci | in what follows since for any given asymmetry parameter £ = 02/ 'o~i, 
we have fi2 = m/ICI- The ratio 



"ID 



<zV/fe (39) 



is then the corresponding electrostatic coupling parameter in ID. (On a more formal level, this definition for the 
coupling parameter may be established by looking at the field-theoretical representation (fT5|) and rescaling the coor- 
dinates i/fii.) Note that in ID the coupling parameter is related to the asymmetry parameter by virtue of 
the electroneutrality condition 

Nqe = -{a-i+a 2 ) , (40) 
where N is the number of counterions of valency q, as 
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Thus, it may be expected that in the limit N — > oo, contrary to the 3D case, the mean-field theory corresponding 
to this system becomes exact! On the other hand, the SC description may be expected to follow simply for N = 1. 
These limiting cases will be discussed further in the forthcoming sections and will be compared with the exact results 
and Monte-Carlo simulations. One should bear in mind that the ID system considered here corresponds to a 3D 
system of mobile charged plates (membranes) confined between two fixed planar charged walls and thus, such a ID 
system with even a single "counterfoil" would have a completely meaningful thermodynamic behavior. 

We note also that both the PB theory (S — > 0) and the SC theory (3 — > oo) for uniformly charged plates are one- 
dimensional theories and should remain valid in 3D as well as in the ID case that we are studying here. However, even 
though the mean-field solution depends only on the transverse coordinate, the field fluctuations about the mean-field 
solution will depend on all coordinates and so it is clear that fluctuations corrections to the mean-field contribution 
depend on the dimensionality of the system: they will be different for a ID system than for a 3D system. 



A. Weak-coupling limit: Poisson-Boltzmann theory 



The expression for the partition function Zn in Eq. (| 1 5[) is up to multiplicative constants given by the functional 
integral 



Z N = / d[<p}exp(-S[<f>]) 



where the action S is given by 

S[<t>} = 



dx 



1 



2q 2 e 2 f3 



— — — KexptUf 
dx ' 



id>—5(z) — id)—S(z — L) 
qe qe 



(42) 



(43) 



In the weak-coupling regime, the leading contribution to the partition function comes from the saddle-point config- 
uration, <po(x), of the action in Eq. (|43[) [32l [33|. where the field <fio(x) = iijjo{x) turns out to be imaginary and is 
proportional to the mean-field electrostatic potential. The saddle-point configuration can be straightforwardly trans- 
lated into a solution of the PB equation and corresponds to an exact asymptotic result in the limit S — > [2CII2H l3~i| . 
The PB equation for the potential, ipo( z ), can be written as [1, [IB] 



with boundary conditions 



d 2 Mx) 
da; 2 



dtpo 
da; 
dipo 
da; 



-q e fine 



-aiflqe = -, 
A* 

R 2C 

<j 2 pqe = 

A* 



(44) 



(45) 



stemming from the electroneutrality of the system. 

Integration of the PB equation gives rise to the first integral of the system of the form 



1 



2q 2 e 2 f3 \ dx 

where the constant po is nothing but the mean-field PB pressure acting between the bounding surfaces Q and 



Po{x) 



k e 



(46) 



(47) 



is the PB number density profile of counterions between the surfaces. As remarked earlier in section ITTll the choice for 
k is arbitrary since it corresponds to a choice of origin for the potential ip(x). In section UIT1 we couched the Fourier 
solution in terms of Fourier coefficients a n = n n b n with the choice n = 1, and so we adopt this choice here, too. 

The nature of the solution ip (x) obviously crucially depends on the sign of the pressure po [10, El[ . Different forms 
are obtained for positive and negative pressures, corresponding to repulsion and attraction between the bounding 
surfaces respectively and were derived previously [42| . 
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Guided by the contact- value theorem Eq. (j27|) , which we emphasize is an exact result, we express the Poisson- 
Boltzmann, or mean-field, pressure in rescaled units. The repulsive PB pressure is given by 



2p /3p 4s . . ,.„, 
Po = — = 7 — rra = a ■ 48 
A Oi/e) 2 



and all the other quantities have been defined above. Furthermore, a — ap and its value is given by the solution of 



tan(aL) = ?f±%. (49) 



For C < 0, and only then, can the PB pressure be attractive. The attractive PB pressure is given by 



where a is now given as a solution of 



2p f3p l B _ 2 

Po = — = t — / To = • (50) 
a{ CTi/e 2 



coth(aX)=-^|. (51) 



The limiting case of zero pressure can be obtained straightforwardly from either limit. The interaction pressure can 
be obviously computed for any value of the asymmetry parameter, but note that any value of £ can be mapped onto 
the interval — 1 < £ < 1 and the pressure only need to be evaluated in that interval of £ values. Note again that the 
PB interaction pressure is the same for a ID as well as for a 3D system. 



B. Weak-coupling limit: fluctuations 



Fluctuations around the mean field depend on the dimensionality of the system and are different for a ID than for 
a 3D system. In order to proceed, one needs to evaluate the a ppr opriate Hessian of the field action in the partition 
function and study its fluctuation spectrum (see Refs. [32l. [33ll43l| for more details). The Hessian of the field about 
the mean-field value is given by 

6 2 S 



5(p(x)<p(y) 



\(p — itp 



= H{x,x') = 



1 



(3e 2 q 2 



dx 2 



-^+(3e 2 q 2 P0 (x) 



S(x — x'), 



(52) 



Po(x) = cxp(— ipo(x)) is the mean-field PB countcrion density. Hence, 

2a 2 



(3e 2 q 2 p (x) = < 



cos 2 a(x — xq) 



2a 2 



Po > 0, 



Po < 0. 



(53) 



, sinh a(x — xq) 

The corresponding correction, J- 2 , to the free energy of the system comes from the functional integral 



where 



A(L) 



s[ip] 



dxdy 



d[ip]exp(-s[ip}), 



4>(0)=x 



?p(x)H(x, x')ip(x')dxdx' , 



so that the corresponding fluctuation part of the free energy can be obtained in the form 

/3.F2 = -In A{L). 



(54) 



(55) 



(56) 



The functional integral Eq. (|54[) can be evaluated exactly, in two different ways. The first method is based on 
the use of the argument principle [36| converting the discrete sum of eigenvalues of the Hessian operator into the 
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logarithm of the secular determinant V\ of the same operator. The second method is based on the Pauli-van Vleck 
approach to calculating the functional integral of a general harmonic kernel. Both results are the same. 

In the first approach the trace-log of the Hessian can be written equivalently in the form that was derived for a 3D 
but can be used in a trivially modified form also for the ID case under consideration. It gives 



case 32, 3 



f3T 2 = hTr \nH(x,x') = |ln 



2l 

2V 



(57) 



where T>\ is the secular determinant, i.e. the determinant of the coefficients corresponding to the appropriate boundary 
condition in the solution of the eigenvalue equation, that can be derived from solutions of 



A/3eVpo(x) /aO) = 



(58) 



One can now write A(L) for the quotient CPi/Pp) -1 / 2 , since the secular determinant depends explicitly on the value 
of the inter-surface spacing, L. Using the fact [32|, [H, [43| that the two linearly independent solutions of Eq. (|58| for 
A = 1 in the repulsive regime are 



f^'(x) = tan ax and f^ 2 '(x) = 1 + ax t an ax 
the secular determinant for the symmetric case C = 1 comes out as 



T>i(L) = a sec 2 



tan a 



a— sec 
2 



and the corresponding free energy contribution from the quadratic fluctuations is thus given by 



1 



— In < a sec a— 



L 



L 

tan I a— 



L 



a— sec a— 



L 



(59) 



(60) 



(61) 



The fluctuation free energy can be regularized so that all irrelevant constants, i.e. all the terms not depending on 
the separation between the bounding surfaces, are dropped, amounting to a rescaling T 2 (L) — > J~i(L) ~ —* oo). 
This corresponds to a subtraction of the part of the free energy for two separate boundaries at infinite separation 
from the total free energy 

In the second method we compute the fluctuations about the mean-field path using the Pauli-van Vleck approach 
[3 EE 13 • As before the thickness of the film is L we take the leftmost and rightmost points of the film to be at 
x = and x = L respectively. Clearly the action of the classical path minimizing s in Eq. (|54|) is a quadratic function 
of the initial and final points of the fluctuating field tp. The generalized Pauli-van Vleck formula tells us that 



i>(L)=y 



ip(0)=x 



d[ip] exp(-s[i/>]) = 



1 ds c [x,y] 
2ir dxdy 



exp(-s c [x,y]) , 



where s c is the classical action minimizing s. As the action is quadratic we may write 

*c[x,y] = \ [f(L)x 2 +g(L)y 2 - 2h(L)xy] . 
In the case of a symmetric charge distribution we have f = g and thus the fluctuation term is given by 

A(L) = 



Solving the equations of motion we find that 
(3q 2 e 2 f(L) = (3q 2 e 2 g(L) = 
pq 2 e 2 h(L) 



2irh{L) V 
P(L)-h*(L)J 



<sec 2 (<xL/2) a[tan(<xL/2) + (aL/2) sec 2 (aL/2)] 



(62) 



(63) 



(64) 



2tan(<xL/2) 
a 



2[1 + {aL/2) tan(aL/2)] 



2tan(aL/2)[l + (aL/2) tan(aL/2)] 
Putting all this together then yields 

A(L) = 



asec 2 (aL/2)[tan(aL/2) + (aL/2) sec 2 (aL/2)} 



(65) 



(66) 
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and the contribution to the free energy due to the fluctuations about mean-field free energy is again obtained, up to 
irrelevant constants disappearing upon rcgularization, as in Eq. (|6ip . 

As we already stated the fluctuation part of the interaction free energy for 3D and ID are completely different (due 
to fluctuations in the plane of the film in 3D which are not present in ID), though the mean- field Poisson-Boltzmann 
result is exactly the same. In particular, the scaling of the fluctuations contribution to the interaction pressure, say 
P2, with the inter-surface distance L turns out to be very different. In 3D, one finds that [421 ] 

p 2 (L) ~ - S 3D (67) 

whereas in ID we find 

P2 (L) ~ - Sid j (68) 

for same-sign surfaces (£ > 0) at sufficiently large separations. Note that the coupling parameter S is defined differently 
in ID, Eq. (|39|) . and 3D, Eq. ((38)) . As one goes to the weak-coupling limit in ID for N — > oo the fluctuation term 
obviously makes a vanishing contribution to the total pressure and is thus not particularly important. The logarithmic 
dependence in 3D is a direct consequence of the contribution from the in-plane modes. In this case, the mean-field 
pressure scales with the inverse square of the separation as Po(L) ~ 1/L 2 irrespective of the dimensionality of the 
system. 

Note that the fluctuations part is always attractive, reflecting the fact that electrostatic correlations mediated by 
counterions always favor attraction between the charged boundaries. Within the WC analysis the fluctuations are 
always assumed to be small as compared with the leading order mean-field contribution. Thus, the total pressure, 
Po + P2, is dominated by the mean-field contribution (except at the equilibrium point where po = (42J) and in 
particular for same-sign surface charges, it still remains repulsive. 

It is also interesting to note that the magnitude of fluctuation pressure relative to the mean-field pressure decreases 
with the separation distance L in 3D, while in ID, it increases with L. For same-sign surfaces, the ratio |p2 /Po | scales 
as ~ InL/L in 3D, while it scales as ~ L in ID. This indicates qualitatively different distance-dependent behaviors 
for the fluctuations and thus, qualitatively different regimes of validity for the loop-expansion approach in ID and 
3D. The latter is determined by assuming that \p 2 /po\ <C 1. Hence, the weak-coupling validity criterion for same-sign 
surfaces in 3D reads 

E W < ±- L . (69) 

This means that at a given coupling parameter, the WC analysis becomes increasingly more accurate at larger 
separations, while as the surfaces get closer a smaller coupling parameter needs to be chosen. In ID, one has 

Sid < p (70) 

which indicates the opposite trend. In particular, for given £ and noting from Eq. (T4"Tj) that Sid = (1 + Q/N, the 
WC analysis becomes increasingly more accurate as the separation decreases for fixed N, or as N increases for fixed 
separation. 



C. Strong-coupling limit 

The strong-coupling approximation coincides with the lowest order of a non-trivial expansion of the partition 
function in terms of the fugacities of the counterions. This expansion may be expressed as a 1/S series expansion 
[20l[2l|. whose leading order term (S — > oo) corresponds to the SC theory. We will not delve into the strong-coupling 
expansion in more detail since it has been exhaustively reviewed in the literature 0, 0, H3, H3 • 

At leading order, the SC free energy is obtained as 

T = W - Nk B T In f e -^ Wl+w ^dV, (71) 
where Wq is electrostatic interaction energy of charged surfaces 



(72) 
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with S representing surface area of each plate, and W\ and Wi are electrostatic interaction energies between a single 
counterion and individual charged surfaces, i.e. 

Wi = -Zpx, W 2 = -^(L- X ). (73) 

Since in the strong-coupling regime the free energy is given via simple quadratures, it is much simpler to evaluate 
it than at the weak-coupling level. Defining the rescaled free energy 



T 



(74) 



we obtain 



| = (1 + C 2 )|-(1 + C) Insinh 



d-C)§' 



(75) 



where L = L/ /i. Differentiating the free energy with respect to the surface-surface distance L we get the corresponding 
pressure acting between the bounding surfaces 



P(L) 



|(1 + C 2 ) + |(1-C 2 ) coth 



(i-O 



L 



(76) 



where, following the discussion in section BV A) we have defined p(L) — 2p(L)/a\. Note that the SC pressure can 
become attractive for both like-charged and oppositely charged surfaces which contrasts with the mean-field theory 
that does not allow attraction between like-charged surfaces. This is because of the strong electrostatic correlations 
mediated by counterions between the charged surfaces for S ^ 1 and has been investigated thoroughly before for 
equally charged surfaces in [20L l2ll ] and for asymmetric surfaces in [42j . 



V. NUMERICAL SIMULATIONS 



We next consider Monte-Carlo simulations of the ID system of counterions treated in the preceding sections using 
both exact as well as approximate (limiting) analytical approaches. Monte-Carlo simulations enable us to access the 
parameter space inaccessible to the limiting WC and SC theories, and thus may be compared directly with the exact 
solution presented in Section 

We proceed by simulating a system of N counterions in a finite interval x G [0, L] in the canonical ensemble by 
applying a standard Metropolis algorithm. The Hamiltonian of the system is defined in Eq. |T]) and the electroneu- 
trality condition is imposed via Eq. (|40|) . We run the simulations by up to 10 s — 10 9 Monte-Carlo steps per particle 
with 10 7 — 10 8 steps used for relaxation purposes. 

The interaction pressure between the two charged boundaries is calculated via the contact- value theorem (see also 
Eq. (27|)) 

fo = P(0) - ^ = P(L) - (77) 

where p(0) and p(L) represent the counterion density at contact with the first (<ti) and the second (172) plate, 
respectively. The resulting pressure computed from the contact density at either surface is found to be the same 
within the numerical errorbars. Simulations were conducted in rescaled units for various number of counterions 
N = 1, . . . , 20 and for the asymmetry parameters £ = —0.5, 0.5 and 1.0. The coupling parameter in each case follows 
from Eq. (jIT]) . 

The simulated pressure is shown in rescaled units in Fig. [T] (symbols) for the rescaled pressure p = (ipi^j {a\j e) 2 
as a function of the rescaled distance, L = L/fi, between the charged boundaries. As seen, the interaction pressure 
decays monotonically with separation distance L for all values of ( and N. Also in all cases, the simulation data 
(symbols) and the exact results (thin solid curves) are nicely bracketed by the mean-field PB result (thick solid curve) 
and the SC result (dashed curve). 

For the symmetric case £ — 1 , the simulation results are spot on the SC curve (dashed line) for N — 1 and approach 
slowly the attractive asymptotic SC pressure at large separations L — > 00, that is p^ — — C 2 - As the number of 
counterions increases the pressure becomes less attractive and eventually for sufficiently large N, the simulation data 
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FIG. 1: (Color online) Rescaled interaction pressure, p, as a function of the rescaled distance, L, between charged boundaries in 
the ID system of counterions. We show the results for three different values of the asymmetry parameter £ = 1.0, 0.5, and —0.5 
(top to bottom). Thick solid lines represent the PB prediction, Eqs. (|48p and (|50p . and the dashed lines are the SC prediction, 
Eq. ([76}. Symbols correspond to MC simulations data and thin solid lines are the exact results, Eq. (|34|l . at different number 
of counterions as indicated in the graphs. 



tend to the PB curve (thick solid curve) exhibiting repulsive pressure at all separations. This is because for large N 
the system effectively splits into two nearly electroneutral halves, where each bounding surface is neutralized by its 
corresponding layer of counterions. This also indicates that the entropic contribution from counterions which favors 
repulsion becomes important at large N. For finite N, the simulation data is described neither by the PB theory nor 
by the SC theory. In this case, an excellent agreement is found between the data and the exact results given by Eq. 
(shown by thin solid lines). 

Similar trends are observed for asymmetric systems as shown for £ = 0.5 and £ = —0.5 in the figure. In the case of 
oppositely charged boundaries (£ < 0), it turns out that the PB and SC curves and hence the exact results roughly 
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coincide and become less distinguishable. A reasonable explanation for this would be that for oppositely charged 
surfaces the counterions mostly feel the effect of the strong uniform external field provided by the surface charges, 
which acts similarly in the strong as well as the weak-coupling limit. Thus the mean-field and the strong-coupling 
approaches should converge. For small external field, as in the case of similarly charged surfaces, the mean-field theory 
depends more on the local counterion density whereas the strongly coupled counterions still feel mostly the external 
field. Thus the difference between the WC and the SC frameworks in the £ < and C > cases. We also emphasize 
that the above discussion holds in rescaled representation as the pressures are plotted here in rescaled units; in actual 
units, Fig. Q] corresponds to different ranges of separation, L, for the WC and SC regimes as the Gouy-Chapman 
length, fi, is typically very different between the two limits, i.e., it would be small at high couplings and large at small 
couplings as may be realized, e.g., by changing the counterion valency at fixed surface charge densities and Bjerrum 
length. 

Note that for oppositely charged boundaries, both the PB and the SC pressure can be attractive at large separations, 
and hence the pressure at any finite N can be attractive. It is also notable that for charged surfaces of opposite sign, 
the PB analysis in general performs much better than for the surfaces of equal sign and that the asymptotic SC value 
is approached more quickly. 

In all cases considered here the theoretical and the simulated values of the interaction pressure converge for very 
small surface-surface separations. In fact, the SC and PB results coincide in the leading order as the rescaled distance, 
L, tends to zero, since both are dominated by the osmotic pressure of counterions. One thus finds 

p(L) ~ -i- L « 1, (78) 

which is of course nothing but the ideal-gas osmotic pressure of counterion confined between the two plates (i.e., 
p = Nk^T/L in actual units). 

VI. CONCLUSIONS 

We have analyzed the statistical physics of an overall electroneutral system, of one-dimensional counterions confined 
between two charged surfaces. This can be seen as a simple model for charged lipid multilayers or stiff mica sheets 
neutralized by counterions and confined between external charged plates. The thermodynamics can be solved exactly 
and a widely varying behavior of effective interaction between the confining plates can be discerned. 

Apart from performing exact statistical mechanics and MC simulations, we have also analyzed these systems within 
the mean-field (Poisson-Boltzmann) or weak-coupling approximation and within the strong-coupling approximation. 
These approximations are in fact independent of the dimensionality of the problem and it is therefore interesting 
to see how they compare with our exact results in one dimension. For a large number of particles it is found that 
the mean-field approximation works well. This is physically understandable as the counterion distribution can be 
reasonably assumed to have a continuous profile as predicted by the mean-field theory and also the inter-particle 
interaction in one dimension is long range. In the case of few counterions it works less well as the effect of correlations 
becomes important. In this regime the SC approximation will perform better as it takes into account properly the 
discreteness of the counterions charges. Indeed, as the strong-coupling expansion is basically a form of the virial 
expansion it is to be expected that it works better in the case of few counterions. 

On the level of the approximate evaluations of the partition function, only the fluctuation contribution to the weak- 
coupling limit depends on the dimensionality of the problem because of fluctuations in the plane of the film in 3D 
that are not present in ID for which the summation over the in-plane modes is absent. In particular, the fluctuation 
contribution to the interaction pressure scales with the distance L as P2 ~ -Sm/I in ID and as P2 ~ — S3D InL/L 3 
in 3D and for same-sign surfaces (£ > 0) at sufficiently large separations. The fluctuational contribution to the 
pressure is always attractive, reflecting the fact that electrostatic correlations mediated by counterions always favor 
attraction between the charged boundaries. Since the coupling parameter in ID, S = Sid, depends inversely on the 
number of counterions, the fluctuation contribution to the interaction pressure that scales linearly with 2 is in this 
case vanishingly small, a situation clearly confirmed by exact and MC evaluation of the partition function. We can 
solve for pi{L) exactly in ID and 3D and evaluate the relative importance of the contributions from fluctuations and 
from mean field theory. For example, in ID with N = 20, C = 1 (Sid = 0.1), we find P2/P0 5; 0.1 for L < 5; whereas 
in 3D we have P2/P0 ^ 0.1 for L > 0.05 when S3D = 0.1, and L > 20 when S3D = 1.0. This confirms the results 
inferred from the asymptotic behaviour of p2(L). 

In the case of symmetric surface charges the mean-field predictions always give a positive (repulsive) pressure 
between the two plates. However the SC and exact results do predict a possible attraction at large inter-plate 
separations. This effect is basically due to strong correlations arising when the counterions cannot neutralize both 
surface charges simultaneously and thus the residual surface charges are both attracted to the residual counterion 
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charge in between the two plates. As the number of counterions increases this effect disappears and the system 
exhibits a more mean-ficld-likc behavior. 

As regards future work it would be interesting to see if the one-dimensional theory can be adapted to model a full 
three dimensional system. There is some hope that a proper reformulation of the statistical mechanics of this system, 
by isolating explicitly the dependence of the local fields in the transverse as opposed to the longitudinal direction 
with respect to the bounding surface normals, might open up new prospects for approximations that would work 
also in the gray area where the WC and the SC fail. This approach may work due to the fact that the mean-field 
and strong-coupling approximations are dimensionality independent. Whether or not it can be expected to work will 
depend on whether charge correlations in the direction perpendicular to the film dominate the physics or whether 
lateral correlations (such as the formation of locally two dimensional Wigner-crystal-like structures in the plane of 
the film) dominate. 
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